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ABSTRACT 

Radiation hydrodynamical transport of ionization fronts in the next generation of cosmological 
reionization simulations holds the promise of predicting UV escape fractions from first principles as 
well as investigating the role of photoionization in feedback processes and structure formation. We 
present a multistep integration scheme for radiative transfer and hydrodynamics for accurate propa- 
gation of I-fronts and ionized flows from a point source in cosmological simulations. The algorithm is 
a photon-conserving method which correctly tracks the position of I-fronts at much lower resolutions 
than non-conservative techniques. The method applies direct hierarchical updates to the ionic species, 
bypassing the need for the costly matrix solutions required by implicit methods while retaining suffi- 
cient accuracy to capture the true evolution of the fronts. We review the physics of ionization fronts in 
power-law density gradients, whose analytical solutions provide excellent validation tests for radiation 
coupling schemes. The advantages and potential drawbacks of direct and implicit schemes are also 
considered, with particular focus on problem timestepping which if not properly implemented can lead 
to morphologically plausible I-front behavior that nonetheless departs from theory. We also examine 
the effect of radiation pressure from very luminous central sources on the evolution of I-fronts and 
flows. 

Subject headings: cosmology: theory — early universe — H II regions: simulation 



1. EARLY REIONIZATION 

The 2003 WMAP discovery of large o ptical depths 
to e le ctron scatteri n g at z ~ 15 ijKogut et al. I 
120031: iSpergel et al. I l2003j) motivated several nu- 
merical and analytical studies of early reionization 
jHaim an fc Ho lder 2003; Ccn 2003: Somcrvillc & Livi^ 
20031 tS okasian. Abel. Hcrna uist. fc SpringeJ 

2003t iCiardi. Ferrara. fc White | 2003t 

Ciardi. Ferrara! MarrL^ Raimondol l2001t 
Wvithe fc Loebl I2003D . The goal of these models 
was to extend the classical picture of reionization to 
account for the large electron fractions at high redshifts 
while still accommodating lower redshift observations. 
In particular, Sokasian, et al. and Ciardi, et al. employed 
large scale cosmological simulations combining gas and 
dark matter dynamics with radiative transfer to follow 
the growth of ionized regions in the early IGM. Global 
electron densities in these calculations were integrated 
over redshift to compute the Thomson scattering optical 
depth Te for a variety of scenarios. 

1.1. Limits to Current Models 

All these simulations assume several free parameters: 
UV escape fractions, Pop HI stellar masses, the pos- 
itive and negative feedback of one generation of UV 
sources upon the next, and the Pop Ill/Pop II rollover 
in mass spectrum and UV production with redshift. 
Furthermore, there was no hydrodynamic response to 
the energy deposited into the gas by the passage of 
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fronts in these models, seriously altering the evolu- 
tion of the fronts themselves. The studies cited above 
postprocessed successive hydrodynamic snapshots of the 
IGM with radiative transfer without energy deposition 
into the gas to evolve primordial H II regions (see 
iMaselli. Ferrara. fc Ciard i (2003) for a summary of the 
latest algorithms applied to cosmological RT). The chal- 
lenge of the next generation of early reionization sim- 
ulations is to capture the radiation hydrodynamics of 
ionization and feedback physics on small scales to de- 
termine the final sizes and distribution of I-fronts in the 
large simulation volumes necessary fo r statistically accu- 
rate structure formation ijBarkana fc Loeb 2004) . 

Upcoming observations able to discriminate between 
early reionization scenarios underscore the need for ab 
initio simulation of the early IGM, in which reioniza- 
tion properly unfolds over many redshifts and gener- 
ations of luminous objects. 21 cm line observations 
in both emission and absorption (by the Square Kilo- 
meter Array) could yield cosmic electron fraction pro- 
files as a function of redshift, beyond current WMAP 
and upcoming Planck measurements that are limited 
to electron column densities. If foreground contami- 
nation can be overcome, these observations might also 
unveil the size and morph ologies of early H II regions 
(jFurlanetto fc Briggjl2004() . Signatures of Pop III stars 
manifest as excesses in the near-IR cosmic background 
may soon be measured by balloon and satellite missions 
(^Santos. Bromm. fc Kamionkowski"2002^ iCoorav et al. I 
2004: Coorav fc Yoshida 2004). JWST wiU also open the 
first direct observational window on protogalaxies with a 
few thousand Pop III stars at 15 < z < 20. 

1.2. Motivation for Cosmological Radiation 
Hydrodynamics 
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The escape of ionizing UV photons from primordial 
minihalos and protogalaxies is mediated by the hydro- 
dynamical transitio ns of ion ization fronts on sub-parsec 
scales (Whalcn, Ab el, fc No rman 2004: Kitavama et al. 
Iin^), and failure to properly capture their breakout can 
alter the final extent of H II regions on kiloparsec scales. 
In some cases static transfer completely fails to predict 
the exit of fronts from early g alaxies by excluding t he gas 
motions that can free them ijWood fc Loebll2000|) . Ra- 
diation hydrodynamical simulations can predict escape 
fractions in the next generation of models by following 
I-fronts as they begin deep within primordial structures 
and blossom outward to ionize the IGM, accurately re- 
solving their true final sizes and morphologies. 

Coupled to reactive networks able to evolve primor- 
dial H2 chemistry, radiation hydrodynamics will also bet- 
ter model the radiative feedback mechanisms known to 
operate in the early universe, which remain to be in- 
corporated in detail in large scal e calculations. Loca l 
entropy injection by UV sources l|Oh & Haimanll2n0,'^ . 
Lyman- Werner dissociation of H2 in minihalos and y^tio- 
togalaxies , and catalysis of molecular hydrogen by free 
electrons lIRicntti. Cnedin. fcSh;^l2nnit fO'Shea et al. I 
120051: iMachacek et al. 112003(1 are key processes governing 
the rise of early star populations and the high-redshift 
ionizing background. On small scales radiation hydro- 
dynamics will also resolve ionized gas outflows in H II 
regions that facilitate the dispersal of metals from the 
first supernovae, exhibit dynamical instabilities poten- 
tially leading to clumping and further star formation, 
and limit the growth of black holes left in minihalos. 
Resolving radiative feedback over a few generations of 
primordial stars will enable their inclusion in large sim- 
ulations over many redshifts with confidence later on. 

On large scales radiation hydrodynamics will be cru- 
cial to determine whether IGM photoheating cascades 
from small to large scales through nonlinear dynamical 
evolution to affect structure formation at later redshifts. 
The expansion of cosmological ionization fronts through 
filaments and voids is also inherently hydrodynamical 
in nature, as is the photocvaporation of minihalos tha t 
can impede these fronts (Shap iro. Iliev. fc Ragall2004|) . 
Static transfer cannot reproduce the outflows confirmed 
by numerical simulations to enhance the photoionization 
of these structures and therefore understimates the ad- 
vance of I-fronts into the early IGM. 

1.3. Overview 

To investigate the numerical issues confronting the 
incorporation of radiation hydrodynamics into future 
large scale structure evolution models we have devel- 
oped an explicit multistep scheme for ionization front 
transport fr om a single po int source in the ZEUS-MP 
hydrocode l|Normanl l2000|) . The issues fall into two 
categories. The first is how to calculate photoioniza- 
tion rates everywhere on the numerical grid, whether by 
ray tracing, variable tensor Eddington factors, flux lim- 
ited diffusion, or Monte C arlo approaches for either sin - 
gle or multiple sources (seeE azoumov fc Ca,rda,11lll2005ll : 
iRiikhorst et al. I l|2005f) for novel raytracing schemes for 
the Enzo and FLASH adaptive mesh refinement (AMR) 
codes). In section 2 we examine the second issue: how 
to couple reaction networks and energy equations driven 
by ionization to the hydrodynamics, which has only re- 



centl y begun to be examined by the cos mology commu- 
nity l)Wehrse. Wickramasinghe. fc Davel l2b0.'7^. Our al- 
gorithm easily extends to multiple frequencies and can 
be readily interfaced with transfer techniques accommo- 
dating many point sources. In sections 4 and 5 we present 
a comprehensive suite of static and hydrodynamic I-front 
test problems utilized to benchmark our code that can 
be applied to validate future methods. In particular, the 
hydrodynamic benchmarks are adopted from an analyt- 
ical study of ionization fronts in power -law density pro- 
files done bv .Franco et al. . (^. 990 ) fsee lYorkd l|1986D for 
a thorough review of numerical and analytical studies of 
classical H II regions). The tests encompass the range of 
I-front dynamics likely to occur in cosmological settings 
and will challenge the versatility and robustness of any 
code. They also expose many features of ionized flows 
that are exhibited by any density gradient. 

We examine in section 6 how individual zones ap- 
proach ionization equilibrium as well as the timescales 
that govern each phase of I-front and ionized flow evolu- 
tion in a variety of density regimes. We discuss how these 
timescales control the timestep advance of the numerical 
solution and explore avenues for future algorithm opti- 
mization. The impact of radiation pressure on I-fronts 
and flows is also reviewed in section 7, and we provide 
an array of improved UV escape fraction calculations for 
Pop III stars in section 8. 

2. NUMERICAL ALGORITHM 

Our modified ZEUS-MP hydrocode solves explicit 
finite-difference approximations to Euler's equations of 
fluid dynamics to gether with a 9-spec ies primordial gas 
reaction network l|Anninos et al. 119971) that utilizes pho- 
toionization rate coefficients computed by a ray-casting 
radiative transfer module. Ionization fronts thus arise 
as an emergent feature of reactive flows and radiative 
transfer in our simulations and are not tracked by com- 
puting equilibria positions along lines of sight, as is done 
in many algorithms. The fluid dynamics equations are 
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-V • (pwjv) - Vp - pV$ - V- Q (1) 
-V • (ev) - pV • V - Q : Vv + Crad 



where p, e, and the v^ are the mass density, internal en- 
ergy density, and velocity of each zone and p (= (7-I) 
e) and Q are the gas pressure and t he von N cumann- 
Richt meyer artificial viscosity tensor l|Stone fc Normal 
|1992() . Crad represents radiative heating and cool- 
ing terms described below. The left-hand side of 
each equation is updated term by term in operator- 
split and directionally-split substeps, with a given sub- 
step incorporating the partial update from the previ- 
ous substep. The gradient (force) terms are computed 
in the ZEUS-MP source routines and the divergence 
terms are calculated in t he ZEUS-MP advcction routines 
UStone fc Norma,niri99^. 

2.1. Reaction Network 
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The primordial species added to ZEUS-MP (H, H+, 
He, He+, He^+, H^, H^, H2, and e~) are assumed to 
share a common velocity field and are evolved by nine 
additional contin uity equations a nd th e nonequilibrium 
rate equations of lAnninos et ai] l)1997j) 

^ = -V . (ftv) + ^ ^ P,,{T)p,p, - ^^P= (2) 

j k j 

where j3jk is the rate coefficient of the reaction between 
species j and k that creates (+) or removes (-) species i, 
and the the ionization rates. The divergence term 

for each species is evaluated in the advection routines, 
while the other terms form a reaction network which is 
solved separately from the source and advection updates. 
To focus our discussion on the gas dynamics of I-fronts 
our present calculations take the primordial gas to be 
hydrogen only, which can be ionic or neutral but not 
molecular. The rate equations reduce to 

= krecnH+TLe - kphUH ~ kcolinHUe (3) 

where k^ec, ^ph, and kco/;, are the rate coefficients for 
recombination, photoionization and electron coUisional 
ionization of hydrogen and the n^ are the number densi- 
ties. 

When the full reaction network is activated we enforce 
charge and baryon conservation at the end of each hy- 
drodynamic cycle with the following constraints 

PbfH = PH + Ph+ + Ph- + PH2+ + PH2 (4) 
1 1 1 

niHTle = Ph+ - Ph- + 2PH2+ + -^PHe+ + l^PHe^ + 
Pb{l- Jh) = PHe + PHe+ + PHe^+ 



where in is the primordial hydrogen fraction, mn is the 
hydrogen mass, and pb is the baryon density evolved in 
the original ZEUS-MP hydrodynamics equations. Any 
error between the species or charge sums and pi, is as- 
signed to the largest of the species to bring them into 
agreement with pi,. 

Microphysical cooling and heating processes are in- 
cluded by an isochoric operator-split update to the en- 
ergy density computed each time the reaction network is 
solved: 

Grad = kphCrUH - Arec»^_f/+ " Apne (5) 
- (Aion + Aexc)nHne 

where kph is the photoionization rate described above, er 
is the fixed energy per ionizati on deposited into the gas 
(set to 2.4 eV as explained in lWhalen. Abel, fc NormanI 
1)200411 ). and Arec, Ac, Aio„, and Aexc are the recom- 
binational, Compton, coUisional ioniza tion, and coUi- 
sional excitation cooling rates taken from lAnninos et al. I 
l)1997ll ). These four processes act together with hydro- 
dynamics (such as adiabatic expansion or shock heating) 
to set the temperature of the gas everywhere in the sim- 
ulations. 



2.2. Radiative Transfer 

The radiative transfer module computes kph by solving 
the static equation of transfer along radial rays outward 
from a single point source centered in a spherical-polar 
coordinate geometry. The fact that the medium usually 
responds to radiation over much longer times than the 
light-crossing times of the problem domain permits us to 
omit the time derivative in the equation of transfer that 
would otherwise restrict the code to unnecessarily short 
timesteps. This approximation is violated very close to 
the central star by rapid ionizations that can lead to 
superluminal I-front velocities. These nonphysical veloc- 
ities are prevented by simply not evaluating fluxes fur- 
ther from the central star tha n light could have traveled 
by th at time in the simulation ijAbel. Norman, fc Madai] 
Il999|) . Our experience has been that this static form of 
the equation of transfer typically becomes valid before 
the I-front reaches the Stromgren radius, and the code 
computes very accurate Stromgren radii and formation 
times. 

Photoionizations in any given cell in principle are due 
to direct photons from the central source through the 
cell's lower face as well as to diffuse recombination pho- 
tons through all its faces. Recombinations within a 
zone occur to either the ground state or to any excited 
state. We adopt the on-the-spot (OTS) approximation 
()Osterbrockll989j) that a case (A) recombination photon 
emitted in a zone is reabsorbed before escaping the zone 
by decomposing both the photoionization and recombi- 
nation terms in the reaction network 

= fc^^'n^+ne -I- fe'-^^n^+rie - feptsTiif (6) 
- kdiffUH - kcoiinnne 

and equating the first and fourth terms, which cancel. 
Here k^"^^ and k^^' are the recombination rate coeffi- 
cients to the ground state and to all excited states while 
kpts and ^diff are the rate coefficients for cell ioniza- 
tions by central stellar photons and diffuse photons. Tak- 
ing ground-state recombination photons to be reabsorbed 
before they can exit the cell guarantees that no photons 
enter the cell from other locations in the problem, reliev- 
ing us of the costly radiative transfer from many lines 
of sight that otherwise would be needed to compute the 
diffuse radiation entering the zone. While the two terms 
are set equal in the equation above, it should be recog- 
nized that ionizations typically occur much more quickly 
than recombinations. No error is introduced in the net- 
work because the much faster photoionizations will sim- 
ply cycle out any ground-state recombinations over the 
timestep the solution is advanced, regardless of the pro- 
cesses' true timescales. 

The OTS approximation is valid anywhere there is a 
sizeable UV photon mean free path across a zone, which 
is the case within the front itself but not usually in the 
hot ionized gas behind it. Case (A) photons emitted from 
an ionized cell might reach the front itself before being 
absorbed because of the very low neutral fraction in their 
path. Such photons do not advance the front, however, 
because the neutral atom they leave behind will remove 
another source or diffuse photon that would have reached 
the front. Since case (A) recombinations globally balance 
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diffuse ionizations they can still be thought to cancel in 
any single zone on average, even in those that are ionized. 

However, there are important scenarios in which the 
OTS approximation fails to reproduce the true ioniza- 
tion of a cosmological structure. For example, dense gas 
clumps in primordial minihalos cast shadows in the UV 
field from the central source. The OTS paradigm would 
produce shadows that remain too sharp by not account- 
ing for the recombination photons that would leak lat- 
erally into these shadows and soften them. Failure to 
capture the correct shadowing may have important con- 
sequences on the growth of ins tabilities expected to de- 
velop in these ionization fronts ijGarcia-Seg^ira fc Francol 
Il996|) . These instabilities are of interest for their poten- 
tial to promote clumping that could later collapse into 
new stars, especially if mixed with metals from a previous 
generation. Methods for accurate and efficient simulation 
of recombination radiation are being studied in connec- 
tion with 3-D simulations of minihalo photoevaporation 
under development. 

Only a single integration of the transfer equation along 
a radial ray from the source at the coordinate center is 
necessary to compute kpts in a cell. The static trans- 
fer equation is recast in flux form for simple solution in 
spherical coordinates: 



V • F 



-XFr 



(7) 
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(12) 



^ rf (0fc+i ~ 0fc) (cosgj - cos6lj-+i) 
hv 

T^ioniz must be converted into the rate coeficient required 
by the chemistry equations. In the ionization term of the 
reaction equation 



fiH = - kph riH 



(13) 



Vph is the ionization rate coefficient, nu is the number 
of ionizations per volume per second, and iiioniz is the 
number of ionizations per zone per second, with the three 
being related by 



# ionizations 
vol ■ sec 



(14) 
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^ioniz is therefore converted to kp/j by 



where x is the inverse mean free path of a UV photon in 
the neutral gas 

X=—_ (8) 



no 



Simple integration yields the radial flux at the inner face 
of each zone on the grid: 



F, 



(9) 



The ionization rate in a zone is calculated in a photon- 
conserving manner to be the number of photons entering 
the zone per second 
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He I and H ionization coefficients as well as Lyman- 
Werner dissociation rates are also evaluated by this pro- 
cedure. 

Our prescription for generating kp/i is photon- 
conserving in that the number of photons emitted along 
any line of sight will always equal the number of pho- 
toionizations in that direction over any time interval. 
This formulation enables the code to accurately advance 
Tfronts with significantly less sensitivity to resolution 
than methods which solve the transfer equation to com- 
pute a zone-centered fiux or intensity to determine kp^: 



_ Fj rf (0fc+i - (/)fc) {cos9j - cosgj+i) 
hv 

minus the number exiting the zone per second 

'^ph, outer — , U 

_ Fi+i r|^i (0fc+i - 0fe) (cos6lj - cos6lj+i) 
hv 

^ Fj e-xt'-'+i-'--) rf i(bk+i - (bk) (cosg^- ^ cos^j+i) 

hu 

Hence, the ionization rate in the cell is 



kph = 4:TTaph— — dv (16) 
hv 

This non-conservative formalism does not guarantee n-, 
to equal iiioniz along a line of sight except in the limit of 
very high resolution. Such methods can require hundreds 
or thousands of radial zones to converge to proper I-front 
evolution, making photon conservation a very desirable 
property ijAbel. Norman, fc Mad au 1999). 

However, photon-conserving methods are not necessar- 
ily resolution independent, as implied at times in the lit- 
erature. For example, if a grid fails to resolve a density 
peak then even a photon-conserving algorithm will over- 
estimate the advance of the front (even in a static prob- 
lem) because of the strong n^ dependence of recombina- 
tion rates on material densities. The resolution needed 
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for correct I-front evolution in hydrodynamical simula- 
tions is governed by the resolution necessary for the gas 
dynamics to converge, since this determines the accu- 
racy of the densities encountered by the front. Neverthe- 
less, photon-conservative schemes are still the method of 
choice in gas dynamical I-front simulations because the 
grid resolution non-conservative methods would require 
to follow the front can be much greater than the resolu- 
tion needed just for hydrodynamic convergence. 

2.3. Network Solution Schemes 
Each equation in the reaction network can be rewritten 



as 



^ Ci{T,nj) - Di{T,nj)n, (17) 

where is the source term representing the formation 
of species i while is the source term describing its 
removal. A fully-explicit finite-differenced solution would 
be 



t+At _ C'At + n, 



1 + D^At 



(18) 



with Ci and (and the Uj comprising them) being eval- 
uated at the current time. A fully- implict update 



t+At _ 



C'+^^At + n,* 
1 + D^+^^At 



(19) 



would take the source terms to be at the advanced time, 
requiring iteration to convergence for the n^ 

The chemistry rate coefficients driving the reaction 
network exhibit a variety of different timescales that 
make the network numerically stiff. The shortest reac- 
tion times in the network would force a purely explicit 
solution into much smaller timesteps than accuracy re- 
quires, while a fully-implicit scheme best suited to the 
solution of stiff networks demands matrix iterations of 
excessive cost in a 3D simulation. The additional costs 
of matrix iteration in implicit schemes can be offset by 
the much longer timesteps they permit because the sta- 
bility of the network is freed from its shortest timescale. 

However, a disadvantage of fully-implicit approaches 
is that they require simultaneous solution of the reac- 
tion network together with the radiative transfer equa- 
tion (needed for the ionization rates in the network) 
and the isochoric energy equation (which sets the tem- 
peratures utilized in the rate coefficients), all eval- 
uated at the advanced time. While such methods 
evolve ionization fronts wi t h very high accura cy in ID 
l|Tenorio-Tagle et al. Ifl98fit iFranco et al. lll99fD . photon 
conservation is sacrificed in the process of achieving 
concurrency between the network, energy, and transfer 
equations at the future time. Hence, these algorithms 
may only achieve their superior accuracy if high prob- 
lem resolutions are employed (8000 radial zones in the 
iFranco et al. I ()199f]j) theory validation runs). Photon- 
conserving fully-implicit stencils which can operate ac- 
curately at much lower resolutions may be possible with 
further investigation. 

Furthermore, while in principle a fully-implicit net- 
work can be evolved over an entire hydrodynamical 
time, this strategy would accrue significant errors in 



many density regimes. iTenorio-Tagle et al. I l)1986j) re- 
ported that in some test cases it was necessary to evolve 
their network by no more than a few photoionization 
timescales in order to compute energy deposition cor- 
rectly. In such instances the additional costs of iteration 
in implicit schemes outweighs their advantage in accu- 
racy over explicit methods because they must ultimately 
both perform a comparable number of cycles to com- 
plete a problem. It should also be noted that highly 
nonlinear and nonmonotonic primordial heating/cooling 
rates have been observed to retard or prevent Newton- 
Raph son convergence in im plicit cosmological calcula- 
tions (|Anninos et al. llT997l) . 

We instead adopt the intermediate strategy of sequen- 
tially computing each n^, building the source and sink 
terms for the i*'* species' update from the i - 1 (and 
earlier) updated species while applying rate coefficients 
evalu ated at the current problem time ijAnninos et al. I 
I1997|) . The order of the updates is H, H+, He, He+, 
He2+, e", H-, H+, and Ha. This approach allows di- 
rect solution of the densities with sufhcient accuracy to 
follow I-fronts in most density regimes with reasonable 
execution times, which are sometimes much shorter than 
for implicit schemes. Anninos et al. found a speedup of 
ten in sequential species updates over an implicit stiff 
solver package in cosmological test cases involving the 
steady buildup of IGM UV fluxes from a metagalactic 
background. 

2.4. Timestep Control 

Two timescales in general govern the evolution of H 
II regions. The many reaction rate timescales can be 
consolidated into a single chemistry timestep defined by 



chem 



He 
he 



(20) 



formulated to ensure that the fastest reaction operat- 
ing at any place or time in the problem determines the 
maximum time by which the reaction network may be 
accurately advanced. The second timescale is the heat- 
ing/cooHng time theat 



the 



^gas 
^rad 



(21) 



connected to the hydrodynamic response of the gas to 
the reactions. The ratio of the two times can depend 
on the evolutionary phase of the H II region or even on 
the current ionization state in a single zone. In gen- 
eral, reaction times are shorter than heating times as the 
ionization front propagates outward from the central UV 
source but cooling times can become shorter than recom- 
bination times after shutdown of the central source. The 
latter circumstance can lead to nonequilibrium cooling in 
H II recombination regions, which can remain much more 
ionized at low temperatures than would be expected for 
a gas in thermodynamic equilibrium, an effect which has 
been observ ed in cosmological H II region simulations 
l|Ricotti. Gn cdin. & ShuU 2001; O'Shca et al. 2005). 

Many strategies have been devised to interleave re- 
action networks and radiative transfer with hydrody- 
namics. After computing the global minimum Courant 
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Fig. 1. — R-type ionization front propagating from a point source in a uniform static medium with no recombinations. Left panel: I-front 
position as a function of time— theory and algorithm for 64, 128, 256, and 512 radial zones. Right panel; neutral fraction profiles-100 Myr, 
1 Gyr, 5 Gyr, 8 Gyr, and 9.7 Gyr. 



timestep for the problem domain, Anninos et al. evolve 
the species in each cell by advancing the rate equations by 
a tenth of the lesser of the chemistry and heating/cooling 
timescales for that cell until a tenth of the cell's heat- 
ing/cooling timescale is covered. At this point the cu- 
mulative energy gained or lost over the chemistry up- 
dates is added to the cell's gas energy in the microphysi- 
cal heating/cooling substep described earlier, but neither 
velocities nor densities are updated. This cycle is then 
repeated over consecutive heating timesteps in the cell 
until the global minimum Courant timestep is covered. 
Cells with the fastest kinetics require the most chem- 
istry subcycles over a heating time: more slowly reacting 
cells covering their heating time with fewer subcycles are 
quiescent during the subsequent cycles required by the 
faster cells. Likewise, kinetics and energy updates in 
cells traversing the global hydrodynamical time in fewer 
heating cycles are halted over the additional subcycles 
the more quickly heating or cooling cells demand. Every 
cell in the grid undergoes the same number of subcycles 
(which continue until the last cell has covered the global 
Courant timestep), but updates in a given cell are sus- 
pended after it has crossed this timestep. New photoion- 
ization rates are calculated every chemistry timestep by a 
call to a radiative transfer module but the other rate coef- 
ficients remain constant over a heating timestep because 
they depend only on temperature; they are updated at 
the beginning of the next heating cycle with the new gas 
temperature. At the end of the hydrodynamical timestep 
full source and advective updates of velocities, energies, 
and conserved total baryonic densities are performed. 

Although sufficient for the slowly rising UV metagalac- 
tic backgrounds in the Anninos et al. calculations this 
subcycling approach does not accurately simulate the 
growth of ionization fronts. Proper front capturing is 
sensitive to velocities that can build up over the heating 
timestep that are not correctly computed by the Anni- 
nos et al. scheme because it does not update velocities 
until many such heating times have passed. The order 
of execution of our algorithm is as follows: first, the ra- 
diative transfer module is called to compute kp/i via eq. 
(|15|) to determine the smallest heating/cooling timestep 



on the grid. The grid minimum of the Courant time is 
then calculated and the smaller of the two is adopted as 
ihydro] in a 1-D Calculation this is 
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Next, from this same set of kph the shortest chemistry 
timestep of the grid is calculated. 



tc 



0.1 



(23) 



The species densities and gas energy in all cells are then 
advanced over this timestep, the transfer module is called 
again to compute a new chemistry timestep, and the 
network and energy updates are performed again. The 
n^ and energy are subcycled over successive chemistry 
timesteps until ihydro has been covered, at which point 
full source and advective updates of velocities, energies, 
and total densities are performed. A new ihydro is then 
determined and the cycle repeats. 

If the chemistry timestep is longer than the global hy- 
drodynamical time the reaction network is only subcy- 
cled once. This more restrictive hierarchy of rate solves 
and hydrodynamical updates is necessary to compute the 
correct velocities in each zone with the passage of the 
front over the wide range of density regimes discussed in 
section 3. Note that cooling times which are shorter than 
recombination times in the problem are easily handled 
because the code will simply cycle the reactions in each 
cell once over the hydrodynamic timestep. Considerable 
experimentation with alternative hierarchies of kinetic 
and hydrodynamical updates and choices of timestep 
control (some involving photoionization times) proved 
them to be less accurate or robust. 

3. IONIZATION FRONT PHYSICS IN POWER-LAW 
DENSITIES 



iFranco et aT\ l|199(]D performed analytical studies of 
1-D ionization fronts from a monochromatic source of 
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Fig. 2. — R-type ionization front propagating from a point source in a uniform static medium with recombinations. Left panel: I-front 
position as a function of time— theory and algorithm for 64, 128, 256, and 512 radial zones. Right panel: temperature profiles for a 512 
radial zone simulation at 1.77 Myr, 5.31 Myr, 21.2 Myr, 63.7 Myr, and 173.4 Myr. 



photons centered in a radial density profile with a flat 
central core followed by an r^" dropoff: 



nnir) 



He 



if r < Tc 
if r > r^. 



(24) 



They considered photon rates that would guarantee that 
the Stromgren radius R5 of the front if the entire medium 
was of uniform density Uc is greater than ic, noting that 
if R5 < ic that the front would evolve as if it were in 
a constant density. Their analysis of I-front propagation 
down the radial density gradients revealed that there is 
a critical exponent 



^crit 



-\ -1 



(25) 



below which the front executes the classic approach to a 
Stromgren sphere of modified radius 



Ruj — Rs 



3 [r^J 



l/(3-2w) ^ 2w/(3-2w) 

Rs 



(26) 



at which point it reverts from R-type to D-type and con- 
tinues along the gradient, building up a dense shocked 
shell before it. Here Rs is the classical Stromgren ra- 
dius the ionization front would have in a uniform density 
medium 
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1/3 



(27) 
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where Ci is the sound speed in the ionized gas. If = 
1.5 the shock and front coincide without any formation 
of a thin neutral shell as in the lu < 1.5 profiles. When 
1.5 < a; < LOcrit the D-type front will revert back to R- 
type and quickly overrun the entire cloud. Since I-fronts 
ultimately transform back to R-type in any cloud with 
density dropoffs steeper than 1.5 this power constitutes 
the critical point for eventual runaway ionization. As 
expected, if w = then R^^ becomes Rs and R(t) exhibits 
the t*/^ expansion of an ionization front in a uniform 
medium l(Ostcrbroc^ ll98?jl . 

Fronts descending gradients steeper than r""'^''" never 
slow to a Stromgren radius or transform to D-type. Re- 
maining R-type, they quickly ionize the entire cloud, 
leaving behind an essentially undisturbed ionized density 
profile because they exit on timescales that are short in 
comparison to any hydrodynamical response of the gas. 
Completely ionized and at much higher pressures, the en- 
tire cloud begins to inflate outward at the sound speed of 
the ionized material. However, the abrupt core density 
dropoff left undisturbed by the rapid exit of the front 
develops a large pressure gradient because of the equa- 
tion of state in the nearly isothermal postfront gas. The 
sharpest pressure gradient is at the ionized core's edge 
at T ^ Tc- This edge expands outward in a pressure wave 
which quickly steepens into a shock that overtakes the 
more slowly moving outer cloud regions. The velocity 
of this shock depends on the initial density dropoff: if 
LOcrit < LU < 3 then 
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If w < 1.5 < LUcrit the front remains D-type throughout 
its lifetime and continues to accumulate mass in its shell, 
expanding as 



If cij = 3 then 
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Fig. 3. — Ionization profiles for a 512 radial zone simulation at 
1.77 Myr, 5.31 Myr, 21.2 Myr, 63.7 Myr, and 173.4 Myr. 




The front will expand forever in the absence of recom- 
binations but will eventually slow to zero velocity as t 
— > oo. We compare this solution to our code results in Fig 
n]for iiH = 10~^ cm"'^, — 10^^ s^"'^, and outer bound- 
ary of 64 kpc for grid resolutions of 64, 128, 256, and 512 
radial zones. The position of the front is defined to be at 
the outermost zone whose neutral fraction has decreased 
to 50%. The time evolution of the front is shown in 
FigH For a given resolution the algorithm is within one 
zone of the exact solution after 2.0 x 10^ yr. As expected 
for a static homogeneous medium, the photon-conserving 
radiative transfer correctly propagates the R-type front 
independently of numerical resolution. The neutral frac- 
tion profiles are extremely sharp because there are no 
recombinations, and they drop essentially to zero behind 
the front. 

Including recombinations yields the well-known result 
for an I-front in a uniform infinite static medium: 



and if w > 3 



R{t) = i?s [1 - exp{-t/Uec)] 



1/3 



(34) 
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2/{S+2-uj) where Rs is the Stromgren radius 



(31) 



where is the initial core radius and 5 is an empirically 
fit function of lo: 



S ~ 0.55 (w - 3) + 2. 



(32) 



The shock has a constant velocity for Ucrit < w < 3, 
weakly accelerates for = 3, and strongly accelerates if 
w > 3. This is in agreement with what would be expected 
for the mass of each cloud. For lu < 3 the cloud mass 
is infinite so a central energy source could not produce 
gas velocities that increase with time. The cloud mass 
becomes finite just above w = 3, the threshold for the 
ionized flow to exhibit a positive acceleration. 

4. ALGORITHM TESTS: STATIC MEDIA 

We present a series of static Tfront tests of increasing 
complexity for initial code validation. In line with the 
non-hydrodynamical nature of these problems, velocity 
updates were suspended but the heating and cooling up- 
dates in eq © were performed. The energy updates were 
necessary to evolve the reaction rates according to tem- 
perature as well as to regulate the timesteps over which 
the reaction network was advanced. The initial gas tem- 
perature in all the static tests was set to 10 K. 

4.1. Spherical I- fronts in Uniform Media 

The simplest test is an R-type ionization front due 
to a monochromatic point source in a uniform infi- 
nite static me dium in which n o recombinations occur 
f\bcl. Norman ~ Madaul \l99^. The radius of the 
spherical front is easily computed by balancing the num- 
ber of emitted photons with the number of atoms in the 
ionized volume: 



Rs 



1/3 



and 



tree = [aB{T)nH] ^ 



(35) 



(36) 



Taking again hh — 10^ cm^'' and — 10 s^ but 
with an outer radius of 10 kpc, we plot the results of our 
algorithm for 64, 128, 256, and 512 radial zones along 
with the analytical solution in Fig |21 The computed 
curves exhibit excellent agreement with theory, with a 
maximum error of 7.5% between the 64-zone solution and 
eq. (|34|) . The code results are again clustered closely to- 
gether because of photon conservation, and after several 
recombination times they converge to a Stromgren radius 
of 7.72 kpc, within 0.3% of the Rs = 7.70 kpc predicted 
by eq. l|?51) . 

The small departure from theory at intermediate times 
evident in Fig [3 arises because eq. (|34|l assumes a con- 
stant aB(T) throughout the evolution of the front. In 
reality, the H II region has a temperature structure that 
changes over time, as shown in Fig [2 At any given time 
the temperature decreases from its maximum near the 
point source to its minimum at the I-front, and this drop 
in temperature with radius can grow to more than 10000 
K at later times. We adopt an average temperature of 
20000 K for as in eq H34|l . The postfront temperatures 
are greatest near the central star because the gas there 
has undergone more cycles of recombination and pho- 
toionization than the gas near the front. Each cycle in- 
crements the gas temperature upward because lower en- 
ergy electrons are preferentially recombined with a net 
deposition of energy into the gas. Successive profiles 
steadily rise in temperature over time for the same rea- 
son. 
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time [yr) time (yr) 

Fig. 4. — Left panel: approach of the I-front to its final Stromgren radius in a static density profile with nc = 10^ cm~'^, rc = 2.1ol6 
cm, constant ob = 2.4773e-13, and an outer boundary of 5.5el6 cm for 25, 50, 100, 200, 500, and 1000 radial zones. Right panel: I-front 
position as a function of time in the density profile with Uc = 3.2 cm"'', Tc = 91.5 pc, constant ag, and an outer boundary of 0.8 kpc 
for 64, 128, 256, 1024, and 2048 radial zones. The uppermost curve is for a 2048-zone temperature-dependent ag simulation. 



The temperature profiles continue to rise well after the 
Stromgren radius is reached as seen in Fig |21 because 
there are no gas motions in which PdV work can be per- 
formed and because coUisional excitation and ionization 
processes are suppressed by the decline of postfront neu- 
tral fractions with time. Neutral fractions fall as ris- 
ing temperatures slow down recombinations. The ris- 
ing gas temperatures in this static problem eventually 
stall when they sufficiently quench recombinations, well 
before other processes such as bremmstrahlung cooling 
arise. Inclusion of recombinations leads to the ionization 
structure of this static H II region in Fig (which in 
part is determined by the temperature profile), in con- 
trast with the very simple neutral fract ion profiles of the 
previous test. iShapiro fc Girous^ l)1987|) extend this clas- 
sical H II region problem to a ionizing point source in a 
uniform medium undergoing cosmological expansion in a 
Friedmann-Robertson- Walker universe, but this test can- 
not be performed by ZEUS-MP at present because scale 
factors are not implemented in the code. 

4.2. I-fronts in Density Gradients 

Franco et al. studied the hydrodynamics of I-fronts in 
power-law gradients but not their time-dependent propa- 
gation in static profiles. Solutions for R-type fronts exit- 
ing fiat central cores into r^'^ gradients do exist for static 
media but in general are quite complicated. 

4.2.1. Bounded I-fronts 

IMellema et al. I l|2005D found that in w = 1 Franco 
et al. static density profiles the radius of the front ad- 
vances according to 



R{t) = i?s <^ 1 + Wo 



exp 



rrt 



Est: 



rcc.corc 



(37) 



where Rs=L/K is the Stromgren radius and Wo(x) 
is the principal br anch of the Lambert W function 
l|Corless et al. llT99l . W(x) is a solution of the algebraic 
equation W{x)e^^^^ = x and must be evaluated numer- 
ically. Here, L = / {A'Kncrc) and K = ncrcCas = 



?'c/ircc,corc, whcrc troc,corc = {ncCas) ^ is thc rccombi- 
nation time in the core and C is the clumping factor. Eq. 
(pTfjl describes the approach of the bounded front to its 
modified Stromgren radius Rg with time. 

W(x) in general is multivalued in the complex plane; 
its principal branch Wq (x) is single- valued over the range 
[-l/e,0], monotonically increasing from -1 to over this 
interval. As the time in eq. H37|) evolves from zero to 
infinity the argument of Wo(x) advances from -1/e to 
0, guaranteeing that the I-front expands from zero to 
Rs in approximately twenty recombination times. Sev- 
eral commercial algebraic software packages can compute 
W(x); we instead uti lized a recursive algorithm by Halley 
llCorless et al. 11799^): 



e^j {wj + 1) 



{wj + 2)(wj e"^^ - x) 
+ 1) 



(38) 



Adopting the first two terms in the series expansion for 
W(x) as an initial guess for each value of x 



W{x) 
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+ 1) 
(39) 



this method typically converged to an error of less than 
10~^ in 2 - 4 iterations per point. 

We show in the left panel of Fig 01 the radius of the 
ionization front on uniform grids of 25, 50, 100, 200, 500, 
and 1000 zones with inner and outer boundaries of 1.047 
X 10"'^^ cm and 5.6 x 10^^ cm, respectively. We set Uc 
= 10*^ cm-3, rc = 2.1 x lO^^ cm, as = 2.4773 x 10"" 
(constant with temperature), and = 5.0 x 10^^ s~^, 
which results in a core recombination time of 0.128 yr and 
a final Stromgren radius of 5.03 x lO^^ cm. ZEUS-MP 
converges to within 4% of eq. (p?7jl by 500 radial zones at 
early times and to within 1% past 1.5 yr (on the scale of 
the graph, the 500 zone curve is essentially identical to 
the 1000 zone plot). The disagreement between theory 
and simulation clearly illustrates that photon-conserving 
schemes can fail to properly evolve ionization fronts if 
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Fig. 5. — Left panel: evolution of the R-type front in the density profile with nc = 10^ cm~^, rc = 2.1 X lO'^^ cm, and outer boundary 
of 5.5 pc for for the central ionizing photon rates shown. Right panel: R-type I-front position in the same r~^ density field extended to an 
outer boundary of 1500 pc with 500 ratioed zones. 



important density features are not resolved, in this case 
the abrupt falloff in density at the edge of the central 
core. Nevertheless, the code still exhibits good agreement 
(> 10%) with the analytical result with relatively few 
(50) zones. 

4.2.2. Unbounded I-fronts 

The evolution of an R-typc front in a static uj = 2 
Franco density field in general involves complex Lambert 
W functions with several branches but reduces to the rel- 
atively simple unbounded solution ijMellema et al. 120051) 

R{t)^r,[l + 2tltrec,coref'^ (40) 

provided that 

"-7 = Tc Uc as (41) 

We plot the position of the ionization front in the right 
panel of Fig ^ on a uniform grid for the resolutions indi- 
cated, with outer radius 0.8 kpc, Uc = 3.2 cm""^, rc = 91.5 
pc, and a constant as — 2.4773 x 10"^^, which yields 
a recombination time in the core of 0.04 Myr and fi^ = 
9.55 X lO^'' s"^ The Mellema et al. I p005,) solutions 
assume a specific temperature for the H II region that 
does not evolve with time, so as was set constant in our 
simulation. Our ZEUS-MP results converge to within 
10% of eq. (gni) by 2000 radial zones (and are within 2 - 
3% over most of the time range) with — 9.65 x 10^*^ 
s~^. An interesting aspect of the analytical solution is 
its sensitivity to the requirement that as is constant, as 
seen in the as (T) plot in the right panel of Fig ^ Re- 
combinations are supressed in a static H II region whose 
temperature varies with radius and increases with time, 
freeing the front to advance more quickly than expected 
from eq. (gUIl- 

R(t) is also strongly divergent from eq. H40|l for photon 
rates above or below eq. (|41|l . as seen in the left panel of 
FiglHl We set n^ = 10*^ cm^^^ ^ 2.1 x 10^^ cm, and as 
= 2.4773 X 10~^^ in a simulation with a 172-zone ratioed 



grid with an outer radius of 5.5 pc. These much higher 
and more compact central densities were chosen for their 
relevance to high-redshift cosmological minihalos and are 
more computationally demanding. The ratioed grid is 
defined by the requirements 



— P ^1 — '^outer 'dinner V j 

where n is the number of radial zones and T:inner was 
chosen to be small in comparison to 5.5 pc but large 
enough to avoid coordinate singularities at the origin. 
We applied a grid ratio (3 = 1.03 to concentrate zones at 
the origin in order to resolve the central core and density 
drop. Our choice of parameters sets — 3.844 x 10^^ 
s~^ and tree, core — 0.128 yr. As shown in Fig [S] the 
ionization front in our calculation stalls for this value of 
n-y but agrees with eq. (|in|l to within 5% if we change fi-y 
to 4.314 X 10*9 s-\ a difference of 12%. If the photon 
rate is increased another 2% the front exits the grid much 
faster than predicted by eq. and if it is decreased 

by 2% the front is halted, implying that the rate set by 
eq. H41(l is the threshold for breakout from an r~^ core 
envelope. 

This result is confirmed by substituting eq. 1)41(1 into 
eq. ((35|l to compute the Rs appearing in eq. H25|l . which 
yields uJcrit = 2. This gradient is just steep enough for 
the central flux of eq. H41() to be unbounded, with the 
position of the front evolving according to the power-law 
eq. ((40|l . We note that while the front escapes for oj 
> uJcrit, it will slow down for a; < 3 because the cloud 
mass is infinite but approach the speed of light for uj > 
3 because of finite cloud mass. The eventual slowing of 
the front in an w = 2 cloud for = 5.0 x 10**^ > fijrrit 
is confirmed at large radii in the right panel of Fig [S] by 
extending the profile used above to an outer boundary of 
1500 pc in a 500 ratioed-zone simulation with (3 = 1.01. 

5. ALGORITHM TESTS: HYDRODYNAMICS 

As in the iFranco et al. I l)1990f) 1-D numerical tests, 
a source of ionizing photons with = 5.0x10*^ s~^ 
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Fig. 6. — Flow profiles of I-fronts and ionized shocks. Column 1; density, velocity, and temperature profiles of an ionization front 
advancing along an tt) = 1.0 density gradient. Output times are 8.64e04 yr, 2.73e05 yr, and 5.71e05 yr. Column 2: same profiles, but for 
an ionized core shock in an a; = 3.0 density gradient. Output times are 1.43e04 yr, 3.23e04 yr, and 5.50e04 yr. Column 3: flow profiles of 
an ionized core shock expanding in an a.) = 5.0 density gradient. Output times are 6123 yr, 9380 yr, and 1.23e04 yr. 



was centered in a flat neutral hydrogen core with = 
2.1x10^^ cm, nc — 10^ cm~'^, and an r""^ dropoff' beyond 
Vc- hv was set to 17.0 eV and the grid was divided into 
8077 uniform radial zones with inner and outer bound- 



aries at 1.047x10"'^'^ cm and 5.50 pc, respectively. Al- 
though not necessary to the benchmarks, for complete- 
ness we employ ed the same interstellar me dium (ISM) 
cooling curve of iDalearno fc McCravl (ITqtI utilized by 
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Fig. 7. — Timescale, temperature, density, and velocity profiles of a D-type front in an oj = 1 density gradient at 3.28e05 yr. 



iFranco et al. I l|199CI() in place of the three last terms of 
eq. jnjl. The photon energy and cooling curve together 
set the postfront temperature and therefore sound speed 
Ci of the ionized gas. Our choice of cooling curve and er 
= 2.4 eV led to postfront temperatures of 17,745 K and 
sound speeds of 15.6 km s^^ in all our runs, in contrast 
to the artificially set 10,000 K temperatures and 11.5 km 
sound speeds of the Franco et al. (1990) runs. We 
applied Ci = 15.6 km s^^ to the analytical expressions 
above for R(t) and rc(t). The initial gas temperature in 
all the hydrodynamical tests was set to 100 K. 

It should be noted that the times appearing in R(t) and 
rc(t) are not taken from when the central source switches 
on. R(t) is the position of the ionization front after the 
initial Stromgren sphere of radius R^, has formed while 
rc(t) is the location of the ionized core shock after one 
sound crossing time across the core 

^ !l ^ 1.351 X 10^" s (43) 

Ci 

so these values must be subtracted from the total prob- 
lem time when comparing the formulae to code results. 
Stromgren sphere formation times t f and radii emerging 
from the run outputs are compiled in Table 1 . In all cases 
the code computed Stromgren radii within a zone width 
of the predicted R„ . 



TABLE 1 
Stromgren Radii R„ and 
Formation Times tj 





Kw (cm) 


tf (sec) 


1.0 


5.970el6 


3.142e09 


1.2 


7.290el6 


8.730e08 


1.4 


1.015el7 


1.035e09 


1.45 


1.148el7 


1.102e09 


1.5 


1.337el7 


1.121e09 



The derivation for R(t) does not account for the hydro- 
dynamical details of the breakout of the shock through 
the I-front so the formula is an approximation which im- 
proves as the front grows beyond R^,. Similarly, rc(t) is 
also an approximation which is increasingly accurate as 
the core shock expands beyond Vc- 

5.1. LO = 1: D-Type Front 

Column 1 of Fig El shows the density, velocity, and 
temperature profiles of a classic D-type ionization front 
expanding in an cj = 1 power law density for the problem 
times listed. ZEUS-MP correctly predicts the formation 
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of a D-type front in this gradient; the confinement of 
the front behind the leading shock is apparent from the 
abrupt temperature drop from 17800 K in the front to 
2000 K in the shock. The front decelerates as predicted 
by eq. (|28|l . from 9 km s~^ to 6.5 km over the 
simulation time. The dense shell thickens as the shock 
accumulates neutral gas, and the rapid decline in post- 
front densities as the problem evolves (especially in com- 
parison to the initial profile) reveals the efficiency with 
which the ionized flow evacuates the cloud core. The 
ionized density profile is flat because any density fluc- 
tuations initially present in the isothermal gas become 
acoustic waves that smooth the variations on timescales 
that are shorter than the expansion time of the front 
(which is subsonic with respect to the 17800 K gas). 

The postfront gas remains isothermal even though cen- 
tral densities fall and the front performs PdV work on the 
shock. This behavior can be understood from the vari- 
ous timescales governing the evolution of the gas energy 
density in the ionized flow 

egas = -pV ■ V + Tphi - A (44) 

We plot photoionizational heating, PdV work, radiative 
cooling, and recombinational timescales defined by 

{Tphi,TpdV,TcoohTrec) — 

(^gas ^gas ^gas \ 
Tp/ii ' |pV • v| ' A ' krecriH+neJ 

along with the hydrodynamical Courant time in Fig [7| 
as a function of radius for the uj = 1 flow. The shortest 
timescales are associated with the most dominant terms 
on the right hand side of eq. H44(l . The heating and cool- 
ing times are identical to within 0.1% behind the D-type 
front, with recombination times that are slightly shorter 
and PdV work times that are three orders of magnitude 
greater. Equal heating and cooling times ensure that 
these processes balance in the energy updates and do 
not change the gas temperature, and the much longer 
PdV timescales demonstrate that the work done by the 
flow behind the front is small compared to the gas en- 
ergy. Even though the ionized densities fall over time 
they still cause enough recombinations for the central 
star to remain energetically coupled to the flow by new 
ionizations and maintain the gas temperatures against 
radiative cooling and shock expansion. 

5.2. UJ = 3: R-type Front and Ionized Shock 

The early evolution of the front in the r~'^ envelope 
is shown in the left panel of Fig |S1 with inner and outer 
problem boundaries of 1.047x10^^ cm and 0.3 pc at res- 
olutions of 250, 500, 1000, and 2000 uniform radial zones. 
The curves are ordered from lower right to upper left in 
increasing resolution. The front exits the core in ~ 0.05 
yr through an essentially undisturbed medium since this 
is much shorter than the dynamical time of the gas. The 
plots converge as the core to envelope transition becomes 
better resolved on the grid. ZEUS-MP reproduces the ex- 
pected slowdown of the R-type front in the central core 
and its rapid acceleration as it exits the density gradient, 
having never made a transition from R-type to D-type. 



Notice that beyond 0.1 pc the slopes of all the curves 
are restricted to the speed of light because the static 
approximation of radiative transfer breaks down in the 
rarified dropoff. As explained in section 6, the mean free 
path of ionizing photons abruptly becomes comparable 
to the size of the grid at that radius and an unrestricted 
static approximation would permit nonzero photoioniza- 
tion rates to suddenly appear all the way to the outer 
boundary. As previously noted, in such circumstances 
the code restricts the position of the front to be 

Rfit) < ctproUem (45) 

This requirement prevents superluminal I-fronts over to- 
tal problem times but not necessarily over successive 
timesteps, as seen in Fig |H1 As the front crosses the 
core boundary R/(t) briefly curves upward with a slope 
greater than the speed of light c before it is abruptly 
limited to the speed of light by eq. (|45|l . Having been 
slowed to less than c in the core the front can cross the 
next few zones faster than the speed of light because be- 
cause the total distance the front has traveled over the 
entire problem time will still be less than ctproUem ■ This 
results in the slight unphysical displacement of the front 
upward by ~ 0.1 pc, which can be seen if one visual- 
izes the true curve to continue up and to the right with 
slope c from the point where the computed slope becomes 
greater than c. This error is negligible in comparison to 
the kpc scales on which the front later expands, and the 
code produces the expected approach of the front veloc- 
ity to the speed of light at later times given that the 
mass of the cloud is finite, as observed earlier. In reality 
no cosmological density field decreases indefinitely so the 
front would eventually slow to a new Strom gren radius 
in the IGM l)Whalen. Abel, fc NormiM]l200|) . 

Columns 2 and 3 of Fig depict the ionized flows de- 
veloping in Lo — 2) and uj = b density fields after the R- 
type ionization front has exited the cloud. In both cases 
the departure of the front from the grid is visible in the 
10000 K gas that extends out to the problem boundary. 
Both clouds are ionized on comparable timescales with 
initially flat profiles before much dynamic response has 
arisen in the gas. However, very distinct flows emerge in 
the two profiles. The uj — i field has a shallower drop 
with smaller pressure gradients that drive a shock that 
accelerates weakly but supersonically with respect to the 
outer cloud. A reverse shock develops that can be seen 
in the tapered peak of the velocity profile and in the den- 
sity maximum that falls increasingly behind the leading 
shock as the flow progresses. The supersonic expansion 
of the core does not permit the central densities to relax 
to constant values as in the D-type front, but they are 
fairly flat out to the reverse shock because the shocked 
region is still dynamically coupled to the inner flow. 

At any given moment the postshock gas temperatures 
are uniform in radius up to the reverse shock, but they 
decrease over time as the flow is driven outward. The 
temperatures are almost flat because the gas heating and 
cooling rates are nearly constant in the fairly uniform 
central densities. The temperatures fall with time be- 
cause the central gas densities are lower than in the to = 
1 flow above. There are fewer of the recombinations that 
enable the central star to sustain the postshock temper- 
atures as in the D-type front, so the flow expands partly 
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Fig. 8. — Early I-front evolution simulations with 250, 500, 1000, and 2000 uniform radial zones. Left panel: I-front position as a function 
of time in an a; = 3 gradient. Right panel: I-front transit through an w = 5 gradient at the same resolutions. 
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Timescale, temperature, density, and velocity profiles of the ionized core shock in an ui 



= 3 density gradient at 4.41e04 yr. 
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Fig. 10. — Timescale, temperature, density, and velocity profiles of an ionized hypersonic shock in an a; = 5 density gradient at l.OleOl 

yr- 




Fig. 11. — Row 1 left to right: density, velocity, and temperature profiles of the ionized core shock in an a; = 2.0 density gradient. 
Output times are 1.61e04 yr, 6.04e04 yr, and 9.06e04 yr. Row 2: density, velocity and temperature profiles of the D-type front in an a; = 
1.42 radial density. Output times are 3.67e04 yr, 9.79e04 yr, and 1.83e05 yr. 
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at the expense of its own internal energy. The timescales 
of eq. (|44|l governing the gas internal energy and the 
Courant times are shown in Fig |51 PdV work done by 
the gas clearly outpaces photoionizational heating out to 
the reverse shock with a net loss of internal energy in 
the ionized gas. The suppression of recombinations is 
evident in the recombination times, which are nearly ten 
times the photoheating times. 

5.3. LO — 5: Hypersonic Ionized Flow 

We show the initial transit of the I-front through an 
r~^ falloff in the right panel of Fig|S|for the same bound- 
aries and grid resolutions as in the uj = 3 plots on the left; 
the curves again progress from lower right to upper left 
as the resolution increases. ZEUS-MP again correctly 
predicts the initial slowing of the front in the core and 
its rapid acceleration down the gradient without ever be- 
coming D-type. There are no qualitative differences be- 
tween the fronts in the two gradients prior to their exit 
from the core because their initial conditions are identi- 
cal: the 0-1 = 5 front also crosses the core in ~ 0.05 yr. 
There is likewise little difference between the two fronts 
after their velocities have been limited to the speed of 
light by the code. The two prominent differences are the 
slower convergence of the curves beyond Tc in the uj = 
5 gradient and their greater superluminal velocities be- 
fore being set to c by the restriction algorithm. Both 
trends are expected: the front will descend the steeper 
gradient more quickly before being limited to the speed 
of light and more grid cells are required to resolve the 
falloff. Again, because the uj = 5 cloud has finite mass, 
the velocity of the front becomes c at later times. 

In this profile a much more abrupt isothermal density 
dropoff remains in the wake of the R-type front. Extreme 
pressure gradients launch the edge of the core outward in 
a strong shock that is essentially a free expansion, as seen 
in the velocity profiles in column 3 of Fig |H| The den- 
sity falls off so sharply that recombinations are quenched 
there and the shock expands adiabatically. The shock ad- 
vances so quickly that it becomes dynamically decoupled 
from the flow behind it, so no reverse shock forms and 
all the flow variables remain stratified with radius. Tem- 
peratures rise to 4 X 10^ K in the shock but drop to 500 
K behind it from adiabatic expansion. The sequence of 
velocity profiles confirms that the core shock accelerates 
over the entire evolution time of the problem, as pre- 
dicted by eq. (|31|l . It becomes hypersonic with speeds 
in excess of 400 km s^^ . The central star is least en- 
ergetically coupled to this ionized flow, with much of it 
being driven by its own internal energy. 

Temperatures in the central ionized gas fall with time 
for the same reasons as in the lu = 3 profiles. The energy 
and Courant times for the lo = 5 fiow appear in Fig 1101 
The hierarchy of timescales is similar to that of the a; = 3 
flow up to the position of the shock, at which point cool- 
ing, heating, and recombination times rise much more 
quickly than their lu — 3 counterparts. The steeper jump 
in timescales is due to the strong suppression of recom- 
binations in the highly stratified densities. 

5.4. UJ = 1.5 and 2.0 

ZEUS-MP accurately reproduces the mildly- 
accelerating ionized core shock expected to form in 
r"'^ radial densities as well as the strong adiabatic 



shock predicted for r^^ gradients, with hydrodynamical 
profiles that are in excellent agreement with earlier 
work perform ed by 1-D implicit lagrangian codes 
ijFranco et al. l fl990). In this section we summarize our 
code results for I-fronts and shocks in r^^-^ and r^^.o 
core envelopes and present their density, velocity, and 
temperature profiles in Fig ^] ZEUS-MP verifies the 
theoretical prediction that a constant-velocity ionized 
shock forms in the uj = 2 gradient after the rapid 
departure of the R-type front, as seen in the upper 
middle panel of Fig ^] This ionized flow does not 
steepen into as strong a shock as in the u; = 3 case, 
which is evident from the lower gas velocities and 
smaller spike in the temperature distribution. Timescale 
analysis indicates that heating and recombination times 
are nearly equal in the postshock gas but not to the 
same degree as in the D-type uj — \ front, so expansion 
occurs partly at the expense of the internal energy 
of the gas. Consequently, the level gas temperatures 
behind the shock decrease slightly over time but not 
to the extent in the w = 3 shock. The evolution of 
the postshock gas is intermediate between that in the 
r^^ and r""^ cases. This fiow regime is relevant to 
primordial minihalos: high resolution simulations of the 
formation of these objects yield spherically-averaged 
baryo n density profiles with 2.0 < uj < 2.5 (Abel et aL\ 
I2002() . The correspondence is not exact because oj can 
vary in radius over this range, which is why both shock 
accelerations and decelerations are observed numerically 
in profiles derived from cosmological initial conditions 
(fA^halcn, Abel, & Norman 2004). 

Our code confirms that the lo — l.b front evolves es- 
sentially as D-type but with no shocked neutral shell: 
the front and the shock are coincident and advance at 
the same velocity. Throughout the evolution of this fiow 
the front is precariously balanced at the edge of breakout 
through the shock and down the gradient. This break- 
through is sensitive to small errors in shock position at 
early times; as a result, ZEUS-MP finds that the I-front 
overruns the shock when uo = 1.45 instead of 1.5 (the 
densities, velocities, and temperatures in the second row 
of Fig^^E^re for an r^^-^^ distribution). However, this 
is a relatively small error, and our algorithm captures 
shock and front positions to within 10% of theory for 
u) — 1.4 and 1.6. Although this discrepancy is still un- 
der investigation, we believe it to be due to small errors 
in energy conservation in our eulerian code. In this re- 
spect lagrangian codes would likely enjoy an advantage 
in proper shock placement because of their conservative 
formalism. 

5.5. Hydrodynamical Convergence 

We demonstrate the convergence of our algorithm to 
eq. (|28() for the I-front position R(j(t) in an = 1 dropoff 
as well as to the ionized shock position eqs. (|29|) . (|30|l . 
and (|31|l for w = 2, 3 and 5 gradients in the first two 
rows of Fig ll2l For each density regime we utilized grids 
with 100, 250, 500, 1000, 2000, 4000, and 8000 uniform 
radial zones. The dashed line in each panel of Fig E| 
is the analytical solution. In Fig 112b. the numerical pre- 
dictions for the position of the w = 1 D-type ionization 
front are ordered from lowest to highest resolution from 
the lower right corner upward to the left toward the ex- 
pected solution. ZEUS-MP converges to within 10% of 
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Fig. 12. — Upper and middle rows: convergence simulations with 100, 250, 500, 1000, 2000, 4000, and 8000 uniform radial zones. Panel 
(a); D-type I-front position as a function of time in an w = 1 radial gradient. Panel (b): time evolution of a weakly accelerating a; = 3 
ionized shock. Panel (c): position of a constant- velocity ionized core shock in an a; = 2 density gradient with respect to time. Panel (d): a 
strongly accelerating hypersonic core shock in an a; = 5 dropoff. Third row: runs with both uniform and ratioed zones compared to theory. 
Panel (e): oj = 1 and 3 runs with 250 uniform zones and 172 ratioed zones. Panel (f): oj = 2 runs with 100 uniform zones and 172 ratioed 
zones. 
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eq. lEHl) by 2000 zones and to within 2 - 3% at 8000 
zones (the 4000 zone curve is indistinguishable from the 
8000 zone solution at this scale). Because eq. does 
not account for the R- to D-type transition of the front, 
the agreement improves the further the front advances 
beyond the Stromgren radius at any resolution. 

Convergence was fastest in the case of the constant ve- 
locity ionized core shock in the w = 2 gradient, as shown 
in Figll2b. The order of solutions is again from lowest to 
highest resolution from the lower right upward to the left 
toward eq. (|29|l . Agreement to within 10% between 10"^ 
yr and 10^ yr and 2 - 3% past 10'' yr was achieved with 
just 500 zones. The trend is somewhat different for the 
weakly accelerating ionized shock in the to = 3 panel in 
Figll2b. The numerical models are all somewhat above 
the analytical curve at later times but the 100 and 250 
zone runs are below it at earlier times. The different 
grids are most easily distinguished by their overshoot at 
late times: with the exception of the 100 zone curve they 
converge downward toward eq. (|3()|l with increasing res- 
olution, with the 8000 zone curve dipping slightly below 
it near the end of the evolution time. With 2000 zones 
the convergence is to within 10% of the expected values 
between 5000 yr and 10 kyr and to within 1.5% by 60 
kyr. 

In the highly supersonic uj = 5 core shock runs shown 
in Fig I12ti. the progression of the numerical curves to- 
ward eq. (|31|) is again from lower right to upper left as 
the resolution increases, but they converge to a solution 
that lies well above the analytical result. Although the 
hydrodynamical profiles for the w = 5 ionized flow in 
Fig m are in excellent qu alitative agreement with earlier 
work ({Franco et al. IIT99Q1 , in this gradient the algorithm 
does not accurately compute the shock placement for rea- 
sons that are currently under investigation. This flow 
regime is rather extreme, with very high shock tempera- 
tures that challenge the energy conservation of our code. 
Such flows are unlikely to develop in the 3-D photoevap- 
oration of the high-redshift minihalos that we plan to 
study. In the 2 < w < 3 gradients of relevance to those 
primordial structures ( Abel et al. 2002), ZEUS-MP is 
quite accurate. 

The convergence trends in the early I-front evolution 
studies in the past three sections suggest that a key factor 
in accurate front placement is proper resolution of the 
central core and its initial dropoff. Ratioed grids with 
higher central and lower outer resolutions can capture 
proper shock placement and front positions with far fewer 
zones than uniform grids, as shown in the bottom row of 
Fig El In Fig 1121 ; we plot ionized core shock radii as 
a function of time in both uj — 1 and 3 density proflles 
for a 250 zone uniform grid and a 172 zone ratioed grid 
with P = 1.03 as defined in section the analytical 
prediction in both panels is again the dashed line. We 
show the corresponding simulations for an w = 2 profile 
in Figll2f but with a 100 zone uniform grid. We achieve 
the same accuracy with 172 ratioed zones as with 2000 
uniform zones in the first two regimes and as with 500 
zones in the lo = 2 gradient (the most convergent of the 
four regimes tested). Equal accuracy with large savings 
in computational resources motivated our use of ratioed 
grids in earlier work, with the only sacrifice being the 
loss of some det ail in the hydrodynamical profile s (e.g. 
compare Fig 3 in lWhalen. Abel, fc NormanI l|2004D to the 



ui = 2 profiles in Fig lllfl . 

6. TIMESTEP EVOLUTION 

The physical timescales governing code timesteps as I- 
fronts and fiows evolve on the grid yield insights about 
the processes that drive the fronts and are the key 
to future algorithm optimization. In general, different 
timescales dominate the transit of the front than the 
rise of ionized fiows behind it. We examine how these 
timescales control the advance of the solution in four 
density gradients: uj = 1.0, 2.0, 3.0, and 5.0. As a 
rule, the photoheating timescales in the cell being ion- 
ized while the front remains on the grid determine the 
global timestep over which the entire solution is updated. 
We first describe how cells approach ionization equilib- 
rium at different radii in these gradients. 

6.1. Ionization Equilibrium of a Single Zone 

The ionization fraction Xg of a zone as a function of 
time varies with distance from the central source, ambi- 
ent density, and type of ionization front. In this section 
we first examine the equilibration of a zone in the fiat 
central core (the same in all four density gradients) and 
then analyze three other zones at radii of 0.25 pc, 1.0 
pc, and 4.5 pc, well into the density falloffs in which 
the distinct characters of the fronts emerge. We plot 
the ionization curves for the four radii in Figs 113b .. c, e, 
and Figll4b for w = 1, 3, 5, and 2 gradients, respectively. 
We tally the corresponding consecutive heating timesteps 
iheat (eq. (|23)) in the four cells as they come to equilib- 
rium in Figs 113b. d, f, and Fig 114b (recall that theat is 
the interval over which full hydrodynamical updates are 
performed). 

6.1.1. Central Zones 

The zone we consider in the central core is the fourth 
from the coordinate origin, but the ionization of the inner 
boundary also merits discussion. Chemical subcycling 
over theat is uniquely extreme in this zone because the 
UV radiation entering its lower face initially encounters 
no electron fraction whatsoever. Since Ug is zero ichem 
is extremely small so far more chemical timesteps are in- 
voked to cover theat than are required by the network 
(several thousand). In practice this only occurs on the 
inner boundary because over the 80 - 90 heating times 
required to ionize the inner zone the hydrodynamic up- 
dates advect a small electron fraction into the next zone, 
substantially reducing the number of initial chemical cy- 
cles there. The network subcycles ~ 100 times in the 
fourth zone in its first heating timestep but only two or 
three times by the fifth timestep. Roughly 100 theat are 
necessary to ionize this zone. 

As shown in Figll3tl. heating times are at first ^ 1 s 
in the central zone but increase to 1 x 10** s by the time 
it reaches Xg = 0.999. The heating rate Crad is great- 
est when the front enters the cell because kp/iU^fer is 
at its maximum but cooling rates are at their minimum 
because of the low initial temperature of the cell. This 
low temperature implies that e^os is also at its minimum, 
ensuring that the first theat is the smallest. As the tem- 
perature and gas energy rise with Xg, cooling processes 
are activated and decrease e-rad, increasing theat- When 
Xg exceeds 0.999 heating times can abruptly dwarf the 
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Fig. 13. — Ionization fraction and t^eat a-s a function of time in zones in the core and at r = 0.25 pc, 1.0 pc, and 4.5 pc as they approach 
ionization equiUbrium. Panels (a), (c), and (e): ionization fractions in lu = 1, 3, and 5, respectively. Panels (b), (d), and (f): the evolution 
of tfteat as a function of time in oj = 1, 3, and 5, respectively. 
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light-crossing time tic of the ceh because coohng sud- 
denly balances photoionizational heating, driving Crad 
sharply downward. The code solution consequently takes 
too large a timestep forward, over which ionization that 
should have commenced in subsequent zones cannot, un- 
physically stalling the young I- front. This only occurs in 
core zones close to the source because tioniz 3> t/c fur- 
ther from the center. Ionization begins in the next zone 
before the previous zone can come to equilibrium, de- 
creasing the code timestep dt and preventing anomalous 
jumps of the solution forward in time. The jumps in cen- 
tral zones are easily remedied by limiting theat to be less 
than a tenth of the light-crossing time there. After a core 
zone is fully ionized the code solution advances by con- 
secutive 0.1 tic steps until eq. H45() permits ionizations in 
the next zone. Hence, the UV flux in the core zones is a 
step function in time, immediately rising to full intensity 
because intervening zones are completely transparent. 

6.1.2. uj ^ 3 and 5 

The heating time profiles for R-type fronts in w = 3 
and 5 gradients appear in Figs ll3t l and f. As noted ear- 
lier, these fronts rapidly accelerate and must be limited 
to the speed of light by the code. The UV flux in the 
outer zones is again a step function in time because, as 
observed above, the front arrives at the next zone before 
the current zone is completely ionized. In this instance 
the step in the radiation field is not quite to its full in- 
tensity because one or more of the intervening zones is 
not fully transparent. As in the core zone, we begin the 
tally of tfieat in the outer three test zones when their 
photoionization rates switch on. The heating times are 
again comparatively small at first and grow by several 
orders of magnitude as the cells become ionized. The 
heating times as a rule increase with distance from the 
source because the smaller outer fluxes generate lower 
photoionization rates. 

The static approximation is clearly violated near 
the inner boundary because the 1 x 10* s ionization 
timescale of the core zone is comparable to its 7 x 10* s 
light-crossing time. The approximation also fails in the 
outer zones of these two gradients even though ionization 
times are much longer than light-crossing times there. 
The explanation lies in the photon mean free paths: as 
shown in Table 2, they exceed the length of the grid as 
the densities plummet with radius. If the widths of the 
fronts overrun the outer boundary the static approxima- 
tion would allow nonzero ionization rates in zones that 
could not have been reached by light by that time in the 
simulation. We must again employ eq. 1)45(1 to prevent 
superluminal velocities in the outer regions. IGM mean 
densities prevent runaway I-fronts in realistic cosmolog- 
ical conditions. In compiling heating time profiles of the 
outer zones we terminate the advance of the front be- 
yond the zone in question to avoid the downshift in code 
timestep associated with ionizations in a new zone before 
the current cell has equilibrated. 

6.1.3. u; = 2 

In Fig 114b the uj = 2 theat profiles at 0.25 pc, 1.0 pc, 
and 4.5 pc have a more complicated structure because 
the UV radiation illuminating those zones is not a step 
function in time. As seen in Table 2 the front is a fraction 



TABLE 2 

UV PHOTON MEAN FREE PATHS. 



radius (pc) 


w = 1 


UJ = 2 


a; = 3 


LU = 5 


0.25 


6.35e-6 pc 


1.29e-4 pc 


4.79e-3 pc 


6.6 pc 


1.0 


1.38e-5 pc 


2.03e-3 pc 


0.299 pc 


5600 pc 


4.5 


6.21e-5 pc 


4.10e-2 pc 


27.14 pc 


1.18c07 pc 



of a zone in width at first but quickly widens to encom- 
pass several zones. Although R-type, the front has a 
velocity well below the speed of light so ionizations pro- 
ceed several zones in advance of the center of the front. 
The radiation field at the leading zone is initially weak 
(having been attenuated by the previous few partially- 
ionized zones) but soon grows to its full intensity as the 
center of the front crosses the cell. As a result the en- 
ergy deposition into the gas is first relatively small but 
increases as the radiation field intensifies. Crad then dips 
back downward when cooling processes are switched on 
as the zone temperature increases. The heating times 
therefore curve downward but then recover upward as 
shown in Fig 114b . The profiles again migrate upward 
with distance from the source. 

The static approximation is well-obeyed in this regime 
given that the ionization times are much longer than the 
zone-crossing times for light and that the speed of the 
front is much smaller than c. The width of the front also 
remains small in comparison to the grid. We begin the 
heating time tally in the three cells when Xg has risen 
to 1 X 10~^ because there is no sudden activation of 
photoionization rates in them. 

6.1.4. Lj = 1 

As we show in Figs II 3h and b, test zones ionized by 
the D-type front in an = 1 gradient exhibit noticeably 
different ionization fractions and heating times because 
the shock reaches them before the front. The shock raises 
their temperatures to a few thousand K and induces col- 
lisional ionizations that leave a residual Xg ~ 0.01 - 02. 
Shock heating lengthens theat when the front reaches the 
cell because Cgas is 200 - 300 times greater than in the 
preshock gas and the initial chemistry subcycling is heav- 
ily reduced by the coUisional electron fraction. The ra- 
diation profile in these cells is a smooth function of time 
because of the UV photons escaping ahead of the front 
into the slightly preionized zone (the width of the front 
itself in these densities is less than a tenth of a zone). 
We initiate the tally of theat in the cell when ionization 
rates surge with the arrival of the center of the front. 

The large fluctuations in heating times at 0.25 pc are 
due to postshock numerical oscillations of the hydro- 
dynamical flow variables. These variations dampen at 
larger radii because the shocked neutral shell thickens 
with time; the densities the front photoevaporates have 
more time to numerically relax to steadier values. Post- 
shock oscillation is also responsible for the early varia- 
tions in the 0.25 pc Xe profile. As in the other gradi- 
ents, theat again increases with distance from the cen- 
tral source. The static approximation is most valid 
in this density regime, with material properties chang- 
ing many orders of magnitude more slowly than light- 
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crossing times. 

6.2. Global Timestep Control: R-type/D-type Fronts 

As the nascent I-front emerges from the central core 
the numerical solution evolves in a succession of heating 
timesteps that rise and then fall as one zone comes to 
equilibrium and ionization commences in the next. As 
discussed earlier, the heating timesteps in the zone being 
ionized rarely achieve their full values in Figs E| and [Bl 
because ionization in the next zone is activated before 
the current zone can come to full equilibrium. Because 
the timesteps are determined by the sum of photoioniza- 
tional heating and cooling in the zone being ionized, they 
are much smaller than Courant times in the postfront 
gas (~ 10 yr) while the front is on the grid. Although 
restrictive, such timesteps are necessary in part because 
of the proximity of the ionized gas temperature to the 
sharp drop in the cooling curve near 10,000 K; longer 
timesteps can cause sudden unphysical cooling in equili- 
brated zones that lose energy in a heating cycle. When- 
ever ionization fronts approach the speed of light the so- 
lution is limited to relatively short steps until the front 
exits the grid (although it should be remembered from 
Figs 1131 and 1141 that the timesteps grow as the radius of 
the front increases). The rise and sharp fall of timesteps 
continues until the I-front exits the outer boundary. In 
oj = 3 and 5 gradients the code expends 85% of its cycles 
propagating the front across the grid (~ 19 yr) and its 
last 15% advancing the core shock to the outer boundary 
(~ 1 X lO'^ yr). 

The transport of R-type fronts in a; = 2 gradients is 
uniquely efficient because of the unusual theat profile of 
the zones. When a cell begins to be ionized the heating 
times are initially fairly long and are again relatively long 
as the zone approaches equilibrium. On average, this 
curve allows the solution to evolve much further in time 
each cycle. 70% of the simulation timesteps are utilized 
for the transit of the front (~ 500 yr) with the remainder 
being spent on the exit of the ionized shock (~ 1.5 x 10^ 
yr). This simulation requires only 20% of the CPU time 
of the Lo = 3 and 5 runs, which require about the same 
time to execute. The much longer heating times in the 
OJ = 1 gradient enable the code to advance the D-type 
front across the grid with only 30% of the cycles needed 
in the a; = 3 and 5 models. However, much longer sim- 
ulation times 6.0 x 10^ yr) are necessary to advance 
the front to the outer boundary because of its relatively 
low velocity behind the shock, so this run executes with 
somewhat greater CPU times than the uj — 2 model. 

6.3. Global Timestep Gontrol: Ionized Flows 

In Figs 114b . d, e, and f we plot global minima of the 
Courant, PdV work, photoionizational heating, recom- 
bination, and cooling timescales along with the code 
timestep dt as a function of time after the fronts have 
exited the grid (except for the ui = 1 case) for w = 1, 
2, 3, and 5 density profiles. After the outermost zone 
is completely ionized the minimum heating and cooling 
times by which the solution is governed settle to nearly 
constant values for the sound crossing time of the cen- 
tral core, approximately 500 yr. The 1 x 10^ - 10^ s code 
timesteps dt are similar to the heating times in core zones 
that have come to complete equilibrium (see Fig I13|) . 
500 yr marks the steepening of the pressure wave at the 



core's edge into the ionized shock that begins to evacuate 
the core. As central densities fall recombinations are su- 
pressed, slowing cooling rates and new ionizations. Pho- 
toheating, recombination, and cooling timescales rise, 
surpassing Courant times at t ~ 1000 yr, at which point 
the code adopts tcour as the new global timestep. The 
noise in the code timestep prior to crossover is due the 
nearly equal heating and cooling rates; their difference is 
sensitive to minor variations in each rate and can modu- 
late relatively large fluctuations in Cgas/e-rad- This noise 
disappears after the algorithm turns to icour to compute 
dt. The minimum iphoto^-recom, and icooi originate in the 
relatively flat densities behind the shock (see FiglHl) and 
continue to rise with the expulsion of baryons from the 
center of the cloud. The three timescales branch away 
from each other at late times in the lo — 3 and 5 pan- 
els because the cooling rates are sensitive to the drop in 
postshock gas temperature (from 18,000 K down to 2000 
K in the lo = 3 gradient) as the ionized flows expand in 
the course of the simulation. The rates do not diverge 
in w = 1 and 2 postshock flow because the gas is either 
isothermal or nearly so. 

The Courant time is approximately 15 yr in the newly 
ionized gas in all four cases but begins to fall after the 
core sound crossing time in the last three regimes as 
a result of the formation of the shock. The minimum 
Courant time coincides with the temperature peak of 
the shock at all times thereafter in the oj — 2, 3, and 
5 flows but resides within the I-front in the lu = 1 gradi- 
ent because the much higher temperature of the ionized 
gas in comparison to the weakly shocked neutral gas. As 
shown in Fig 114b . the Courant proflle is almost flat in 
this flow regime because the postfront gas is isothermal. 
The ionized core shock in the w = 2 gradient has a nearly 
constant velocity after formation so its Courant time re- 
mains flat after falling to 10 yr in Figll4bl. The oj — 3 
flow only weakly accelerates and therefore exhibits a sim- 
ilar profile in Fig 114b with flnal Courant times of 4.5 yr 
in the stronger shock. Courant times in the hypersonic lo 
= 5 flow continue to fall in Figll4F because of its strong 
acceleration, eventually dropping to 0.6 yr by the end of 
the simulation. 

The flrst dip in the PdV timescale at t ~ 480 yr in the 
a; = 2, 3, and 5 plots results is due to outflow of heated 
gas from the inner boundary zone. The oscillations af- 
ter the dip at t ~ 800 yr are associated with the launch 
of the core shock. The minimum work times later track 
the forward edge of the shock at all times thereafter in 
all four regimes. The behavior of the w = 2, 3, and 5 
PdV profiles after the rise of the shock mirrors that of 
the Courant profiles. The constant speed uj = 2 shock 
uniformly accelerates parcels of upstream gas through- 
out its evolution so its PdV timescales remain fairly flat. 
In contrast, the strong w = 5 shock accelerates upstream 
fluid elements to increasingly higher speeds as it advances 
so its work timescales continue to fall. The profile of the 
gently accelerating lu = 3 shock falls in between these two 
cases, sloping downward and then evening out at later 
times. PdV timescales in the D-type front actually rise 
over time because of the deceleration of the shock. The 
momentary spike in the profile at t 2000 yr (also mani- 
fest in the Courant profile) occurs at the inner boundary 
and is likely due to pressure fluctuations associated with 
acoustic waves there. 
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iQ' Itf ifl* ia= ia* 10* 



Fig. 14. — Row 1: same as in FigQSlbut for an = 2 density gradient. Panel (a): ionization fraction. Panel (b): heating/cooling time 
theat- Rows 2 and 3: global minima of the Courant, PdV work, photoionizational heating, recombination, and cooling timescales along 
with code timestep dt as a function of time after the fronts have exited the grids (except in the ui = I case). Panels (c) and (d): w = 1 
and 2 density gradients, respectively. Panels (e) and (f): u! = 3 and 5 density profiles, respectively. 
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Several features in the uj — I plot deserve attention. 
The PdV timescale oscillations at t < 400 yr are from 
sound waves in the flat core. The postshock numeri- 
cal oscillations of the flow variables responsible for the 
rapid fluctuations in the heating profiles continue to 
cause rapid variations in the heating, cooling, and work 
timescales throughout the simulation, dominating the 
code timestep until the average theat rises above the 
Courant time at t ^ 3 x fO^ yr. It is important to 
remember that the noise in theat and tcooi in this panel 
is not readily apparent in Fig |7| because the global min- 
ima in Fig occur at the rear edge of the shock. The 
two timescales quickly become equal behind the shock. 

7. THE ROLE OF RADIATION PRESSURE IN I-FRONT 
TRANSPORT 

Direct radiation forces in astrophysical flows can 
strongly influence their evolution in scenarios like ac- 
cretion onto compact objects or line-driven winds from 
massive stars. In order to assess the radiation pressures 
exerted by Pop III stars upon their parent halos we added 
an operator-split update term to the ZEUS-MP momen- 
tum equation: 

iiP^) = - f x{i^)F{'^)d^ (46) 

We utilized a flux corresponding to the blackbody spec- 
trum of a 200 M0 star with Te// — 10^ K normalized 
to yield the emission rate of all photons above the h^ 
drogen ionization threshold given in Table 4 of iSchaeren 
()2002j) and centered in the LCDM spherically-averaged 
halo discussed in section|Sl The neutral hydrogen compo- 
nent of the extinction coefficient exits the integral, which, 
along with the integral required to normalize the flux, 
only needs to be evaluated once. X varies over ionization 
timescales because of its dependence on the neutral hy- 
drogen fraction, so radiation momentum transfer into the 
gas must be accrued over consecutive chemical timesteps. 
This is done by multiplying the force integral above in ev- 
ery problem zone by the current chemical timestep and 
summing this product over subsequent chemical times 
until the hydro timestep has been crossed. The sum is 
then divided by the total density to obtain the update to 
the gas velocity in each zone. This method neglects spec- 
tral hardening across the front that could only be cap- 
tured by a full multifrequency treatment of the radiation 
transfer, but these effects are likely to be unimportant in 
the high densities of the primordial halo. 

Gas densities and UV fluxes are largest at the center of 
the halo so radiation forces exert their greatest influence 
in the earliest stages of primordial H II region evolution. 
UV photons can accelerate both the small neutral frac- 
tions in the postfront gas as well as the thin semi-neutral 
discontinuity at the front. Gravitational, thermal pres- 
sure gradient, and radiation forces in the early I-front are 
compared at 22 yr and 220 yr in Fig^] Thermal pres- 
sure gradients cancel gravity forces upstream of the shock 
leading the front (because of the hydrostatic initial con- 
ditions of the gas) but dwarf them, as expected, within 
the H II region itself. The radiation force is also much 
greater than gravity within the H II region, an interest- 
ing result given that neutral fractions are less than 10^^ 
there. Thermal pressure dominates radiation in the gas 



on either side of the I-front except in the partially-neutral 
layer of the front where the radiation force spikes be- 
cause the product of the flux and neutral fraction peaks 
there. In contrast, the radiation pressure is somewhat 
larger than thermal pressure in the postfront gas at 220 
yr in this simulation, but their ratio in general is sensi- 
tive to changes in the very small neutral fraction in the 
ionized gas that can occur problem resolution or initial 
conditions. Comparison of the direct UV forces upon 
the postfront gas at 22 yr and 220 yr clearly shows that 
they diminish with time. As the 10^ K gas expands, it 
is driven outward with an accompanying drop in central 
densities that suppresses recombination rates. Central 
neutral fractions fall (eventually to 10~^), with a corre- 
sponding loss of radiation pressure. O ur radiation force 
profiles are in general agreement with iKitavama et al. I 
( 2004), who considered the same 200 Mq blackbody spec- 
trum but applied a different initial density profile. 

Our simulations reveal that the radiation forces in the 
ionized gas behind the front enhance its velocity by less 
than one percent, even though they are comparable to 
the thermal forces there. This is due to the fact that the 
thermal pressure gradients in the shocked gas and just 
behind the front are primarily responsible for the accel- 
eration of the flow, and these gradient forces are large 
in comparison to the radiation forces except across the 
few photon mean free paths (mfp) of the front. Direct 
calculation of the UV acceleration of fluid elements in 
the I-front layer itself is currently impractical for two 
reasons. The mfp of UV photons through the baryonic 
densities typical of 10^ Mq halos is approximately 10^^ 
pc, well below the resolution limit of an eulerian calcula- 
tion that must accommodate the much larger dynamical 
scales of the H II region. Failure to resolve the front can 
cause a code to interpret the very large radiation force 
peak to act upon the entire problem cell over the time 
required to ionize all of the cell, when in reality the force 
spike only operates on the fluid parcels in the extremely 
thin front for just the time required to ionize this layer 
(such an error led to unphysically large gas velocities in 
early trials of our code). Furthermore, the intense UV 
flux in such proximity to the central star can violate the 
static approximation to the transfer equation, necessi- 
tating the use of the fully time-dependent equation to 
ensure accurate transport of the I-front, even across an 
appropriately rezoned grid. 

A simple estimate of the radiative acceleration of gas 
elements in the front is possible by recognizing that it 
acts on the layer for at most a few ionization timescales, 
between 10''' s and 10^ s at the times indicated in Fig 
1151 The product of the radiation force peak in Fig 1X^1 
and these photoionization times yields an upper limit to 
the velocity imparted to the front, which we find from 
the data above to be less than 2 km s~^ , too little to 
significantly alter the evolutionary outcome of the ion- 
ized flow. This estimate is an upper bound because the 
force actually decreases as the parcels are ionized, while 
this procedure takes it to be constant (and at its maxi- 
mum strength). When the thin layer is ionized, the large 
radiation force upon it evaporates and continues on to 
the next layer. The global radiation momentum transfer 
from the thin front into the gas accelerates it by at most 
1 - 2 km s^i . 
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Fig. 15. — Thermal pressure gradient, radiation pressure, and 
gravitational accelerations taken at 22 yr (profiles on the left) and 
220 yr (profiles to the right). 



8. IMPROVED POP III UV ESCAPE FRACTION 
ESTIMATES 

The breakout of UV radiation from Pop III 
minihalos was explored in an earlier paper 
ijWhalen. Abel, fc Normal 1200^ (hereafter WAN04); 
we now extend this study to lower stellar masses with 
improved problem setup. We adopted the spherically- 
averaged LCDM baryon field of a 5 x 10^ Mq dark 
matter minihal o formed at z = 18 f or our initial number 
density profile llO'Shea aL II2005I) in place of the CDM 
density profile o f I Abel p.t al. lEfim . Also, in our earlier 
work, the inner boundary was set to 0.06 pc in the CDM 
protostellar density, regardless of the mass assumed for 
the central star later forming at its center. Launching 
the I-front from this radius would neglect the 100 - 200 
Mq gas mass interior to this radius not consolidated into 
the protostar in the lower stellar mass cases, as seen in 
Fig 1161 The inner boundary was therefore set to be the 
radius enclosing a gas mass equal to the mass of the star 
whose UV escape fraction was being studied. Although 
this inner few hundred Mq is small in comparison to the 
outer tens of thousands of solar masses later ionized in 
our previous studies, they were included because their 
significantly higher number densities and recombination 
rates might disproportionately retard the advance of the 
front. Apart from these two modifications, the problem 
setup was the same as in the WAN04 simulations. 

An important consequence of shifting the inner prob- 
lem boundary in accord with the protostellar mass is 
that the nascent I-front can now encounter much higher 
central densities, up to 10* cm~'^ for the 25 Mq runs 
compared to the lO'' cm~'^ densities in WAN04. The 
accompanying rise in photoionization rates restricts the 
timestep control to much shorter advances, making 25 
Mq the lower practical computational limit in these den- 
sity profiles. Despite gathering interest in the physi- 
cal processes that cut off accretion onto Pop III proto- 
stars (O mukai & Inutsuka 2002; Omukai & Palla 2003; 
iRipamonti fc Abe]ll2004D . they are not understood well 
enough to provide firm estimates on the final masses 
of these objects, only that they are likely 30 - 300 Mq 



TABLE 3 

Final I-front radii, breakout times, and Fesc- 



M. (M0) 


radius (pc) 


^brkout (yr) 


''^brkout (pc) 




25 


5.0e-3 









40 


8.8e-3 









80 


1978 


3.3e05 


7.9 


89% 


120 


2633 


2.3e05 


7.2 


91% 


140 


3186 


1.7e05 


4.9 


93% 


200 


3505 


1.3e05 


3.7 


94% 


260 


3856 


9.0e04 


2.8 


96% 


300 


4441 


8.5e04 


2.6 


96% 


400 


4665 


5.9e04 


2.3 


97% 


500 


5132 


1.9c04 


1.0 


99% 



ijAbel et al. ll2002HBromm et al. Ill999() . Future models 
will also be needed to determine the central envelope den- 
sities and infall through which primordial I-fronts truly 
emerge; the exit of ionization fronts through accretion 
infall will be studied in forthcoming 3-D simulations. 

Results for stellar masses ranging from 25 Mq to 500 
Mq using t he tim e-averaged photon rates from Table 4 
of jSchaererj (|2Q03) 

appear in Table 3. Since no UV pho- 
tons exit the virial radius of the halo before the front 
transforms from D-type back to R-type, fesc is simply 

J. ttnsl threakout (47) 

where tmsi and tbreakout are the stellar main sequence 
lifetime and time to I-front breakout, respectively. Com- 
parison of these final I-front radii with those of WAN04 
for M > 100 Mq indicates that they differ by 20% at 120 
Mq but only by 4% at 500 Mq in a trend toward greater 
agreement with mass. WAN04 excluded the most inter- 
vening mass at 120 Mq and the least at 500 Mq, so it is 
not surprising the I-front is delayed the most in the first 
case. As expected, breakout radii and times decrease 
with increasing stellar mass 

The I-front breakout from the 80 Mq case shows that 
lower-mass Pop III stars still effectively ionize their local 
environments and can establish ionized outflows capable 
of dispersing any elements ejected by pulsational insta- 
bilities before the st ar collapses directly into a black hole 
l)Heger et al. ll2?)03l) . In contrast, the primordial enve- 
lope traps the I-fronts of the 25 Mq and 40 Mq stars 
to form ultracompact H II (UC H II) regions with life- 
times of 6.5 Myr and 3.9 Myr, respectively. The higher 
central recombination rates overwhelm the lower ioniz- 
ing photon rates of a lower mass star to stall the front 
at subparsec radii. A star in this mass range cannot re- 
main on the main sequence long enough for the ionized 
pressure bubble to expand and free the front. 

9. CONCLUSION 

We have presented an explicit multi-step scheme for 
integrating the coupled equations of radiation transport, 
ionization kinetics, radiative heating and cooling, and hy- 
drodynamics. We have validated our method against a 
large battery of analytic test problems for both static and 
moving media. While these are ID tests, the method is 
easily extended to 3D. An important strength of our mul- 
tistepping algorithm is its applicability to a variety of ra- 
diation transport schemes for computing photoionization 
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rates on a grid. Many problems can be addressed by this 
method in its current state as a single source code, with 
more soon being possible with the planned activation of 
the full 9-species network as well as multifrequency and 
multisource upgrades. Full 3-D radiation hydrodynami- 
cal trials of the code are currently underway to simulate 
the escape of UV radiation from high-redshift minihalos 
formed from realistic cosmological initial conditions. One 
important question to be answered is whether the UV 
escape fraction cutoff observed by Kitayama, et al. with 
low stellar masses persists if three-dimensional instabili- 
ties arise at the ionization front. Instabilities may open 
channels out of the halo and permit the exit of UV ra- 
diation that would not otherwise have escaped, enhanc- 
ing UV escape fractions from low mass stars. On the 
other hand, recent high dynamical range AMR simula- 
tions of primordial star formation (O'Shea et al. 2005) 
have revealed the formation of disks around protostars. 
These disks may reduce escape fractions even from high 
mass stars by blocking UV photons along lines of sight 
in the plane of the disk with large optical depths. The 
radiation hydrodynamical evolution of these objects will 
enable the simulation of more realistic initial conditions 
for the explosion of supernovae and subsequent chemical 
enrichment of the early universe. Such simulations would 
also set the stage for self-consistent modeling of the ener- 
getics of cosmological miniquasars an d the propagation 
of their radiation into the early IGM l)Kuhlen fc Madaul 
120051 

While we have modified our algorithm to study the 
photoevaporation of minihalos in the vicinity of a mas- 
sive Population III (P op III) star in 1-D at high redshifts 
ijO'Shea et al. \^009L the ionization of these objects in 
3-D with primordial chemistry remains to be done. The 



catalysis of molecular hydrogen within these structures 
in the presence of SUV and x-ray backgrounds as they are 
being ionized by external fronts is a key issue in radiative 
feedback processes in the early universe that is not well 
understood. Another process within the realm of study 
of our single source algorithm is the potential cutoff of 
accretion onto Pop III protostars by nascent I-fronts as 
the star enters the main sequence. Radiation hydrody- 
namical simulation of accretion cutoff processes with im- 
proved stellar evolution models at their foundation will 
provide firmer estimates of the true mass spectrum of 
Pop III stars, a key ingedient in large scale calculations 
of early reionization. Studies of the escape of UV radia- 
tion from high redshift protogalaxies of ^ 1000 Pop III 
stars are currently being planned in connection with the 
addition of cartesian multisource VTEF radiative trans- 
fer (Paschos & Norman 2005-) to ZEUS-MP. Having ex- 
haustively validated our coupling scheme for radiative 
transfer and hydrodynamics in 1-D with the array of 
static and hydrodynamic tests presented in this paper, 
we can now apply the algorithm to these relevant and 
timely problems in three dimensions with confidence. 
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